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1.0  Introduction 


The  Laser  Doppler  Velocimeter  (LDV)  has  been  a  valuable  tool  for 
experimental  investigation  in  fluid  mechanics.  The  obvious  advantage  of  LDV 
is  its  high  frequency  response  and  its  ability  to  measure  direction  and  magnitude 
of  the  velocity  accurately  and  non-intrusivcly  under  conditions  where  other 
instruments  provide  questionable  results  or  cannot  be  used  at  all. 

Although  commercial  pointwise  LDV  systems  are  available,  further 
instrument  development  is  needed.  In  order  to  obtain  more  detailed  features  in 
certain  complicated  flows,  such  as  turbulent  separated  flows,  "instantaneous" 
velocity  profiles  should  be  obtained.  Therefore  scanning  LDV  systems  are  needed 
because  they  yield  nearly  instantaneous  velocity  profiles  and,  concurrently, 
reduce  the  data  acquisition  time  of  velocity  information  in  at  least  one  direction. 
Furthermore,  three-dimensional  measurements  should  be  made  so  that  important 
space-time  flow-structural  information  can  be  obtained. 
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This  work  is  concerned  with  the  development  of  a  three-component  rapidly 
scanning  laser  Doppler  velocimeter  and  the  associated  data  acquisition. 

1.1  Background 

Laser  Doppler  velocimetry  (LDV)  is  an  optical  technique  for  measuring 
the  local,  instantaneous  velocity  of  a  flow.  The  basic  principle  of  LDV  is  to  split 
a  laser  beam  into  two  equal  intensity  beams  which  are  made  to  intersect  at  the 
measurement  volume.  The  intersecting  beams  form  stationary  interference  fringe 
patterns  inside  the  measurement  volume  formed  by  the  crossing  beams. 

The  fringe  pattern  is  made  of  light  and  dark  regions  and  if  the  laser  beams 
cross  at  their  waists,  parallel  fringes  will  be  formed  with  equal  fringe  spacing. 
The  spacing  of  the  fringes,  ?.AB,  is  given  by 

; 

AB  2  sin  cp 

where  <p  is  the  half  angle  between  the  intersecting  pair  of  beams  and  ?.  is  the 
wavelength  of  the  laser  light. 

In  order  to  measure  the  velocity  of  the  flow,  tracer  particles  must  be 
present  in  the  fluid.  These  particles  can  be  naturally  occurring  or  artificially 
induced  into  the  fluid  but  should  be  of  appropriate  size,  density  and  shape  so  that 
the  particles  follow  the  flow.  When  a  particle  passes  through  the  measuring 
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volume  it  scatters  light  that  can  be  detected  by  a  photodetcctor.  The  Doppler 
frequency  of  the  scattered  light  is  directly  proportional  to  the  particle  velocity 
which  is  given  by 


fd  = 


2\U\ 


sin  (p 


where  |  U\  is  the  magnitude  of  the  instantaneous  velocity  component 
perpendicular  to  the  optical  axis. 

For  more  information  on  the  principles  and  practices  of  LDV  see  Durst, 
Melling,  and  Whitelaw  (1981). 


1.2  Previous  Work:  Scanning  LDV 


A  brief  description  of  a  few  scanning  LDV  systems  developed  and  used  so 
far,  is  given  by  Antoine  (1985)  and  Simpson  (1989).  Bcndick  (1971)  described 
an  on-axis  scanning  LDV  that  used  a  translational  oscillating  mirror.  This 
system  was  used  for  instantaneous  velocity  measurements  in  stead"  and  pulsating 
water  flow  in  a  glass  tube  of  6  mm  I.  D.  The  operation  of  this  design  is  limited 
to  scan  speeds  of  0.4  m/s  due  to  the  inertia  of  the  moving  optics. 

A  two-color  dual-beam  backscatter  LDV  system  accomplishing  a  scan  by 
translating  a  lens  in  the  direction  of  the  optical  axis  was  reported  by  Grant  and 
Orloff  (1973).  Scan  rates  were  limited  to  1.5  m  s  because  of  inertial 
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considerations.  More  information  concerning  the  application  of  this  design  is 
given  by  Orloff  and  Biggers  (1974)  and  Orloff,  et  al.  (1975). 

A  backscatter  scanning  system  was  reported  by  Rhodes  (1976).  It  is  able 
to  scan  a  distance  of  30  cm  at  a  frequency  of  30  Hz,  and  measure  velocities  at  16 
discrete  positions  using  a  large  rotating  wheel  containing  16  ports.  For  more 
information  concerning  this  design,  see  Gartrell  and  Jordan  (1977)  and  Meyers 
(1979). 

An  optical  system  capable  of  measuring  true  instantaneous  velocity 
profiles  was  reported  by  Nakatani,  et  al.  (1978).  Instead  of  using  a  moving 
scanning  device,  it  employed  a  cylindrical  lens  to  form  a  vertical  measurement 
volume  along  a  straight  line.  The  design  is  relatively  impractical  and  expensive 
because  a  scries  of  optical  fibers  connected  to  photodetectors  is  needed  to  collect 
data  over  a  large  scan  range  with  good  resolution. 

The  design  of  Durst,  Lehmann,  and  Tropea  (1981)  used  a  relatively  large 
rotationally  oscillating  mirror  in  front  of  a  conventional  LDV  optical  system  to 
scan  the  measurement  volume  perpendicular  to  the  optical  axis  along  an  arc. 
Mean  and  RMS  velocity  profiles  agreed  well  with  pointwise  measurements  for 
low  scan  frequencies.  The  inertia  of  the  oscillating  mirror  limited  the  scan 
frequency  to  about  15  Hz. 

Owen  (1984)  developed  a  single  velocity  scanning  LDV  and  made 
measurements  in  both  water  and  air  flows.  A  six-sided  rotating  polygon  mirror 
was  used  and  the  scan  rate  was  limited  to  125  Hz  due  to  restrictions  on  the  data 
acquisition  rates  and  the  number  of  points  required  for  one  profile. 
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Improvements  were  suggested  for  a  second  phase  of  the  work  which  will  enable 
two-component  real-time  scans  through  high  speed  air  flows.  The  disadvantage 
of  a  rotating  mirror  polygon  is  that  the  reflected  beams  are  in  unwanted 
directions  90%  of  the  time. 

Chehroudi  and  Simpson  (1983)  developed  a  single-component  rapidly 
scanning  LDV  that  scanned  up  to  150  Hz  over  40  cm.  The  proper  operation  and 
usefulness  of  this  design  was  demonstrated  in  a  separated  flow  studied  in  a 
boundary  layer  wind  tunnel.  "Instantaneous"  velocity  profiles  could  not  be 
obtained  due  to  the  data  acquisition  system  used. 

Antoine  and  Simpson  (1986)  extended  the  design  above  to  a 
three-component  system.  A  Ronchi  ruling  was  used  to  produces  fringe  patterns 
in  the  horizontal  plane  to  measure  V.  The  fringe  patterns  produced  by  the 
Ronchi  ruling  were  unsatisfactory. 

Econonou  extended  the  design  of  Chehroudi  and  Simpson  (1983)  to  permit 
scans  in  the  plane  perpendicular  to  the  mainstream  flow  velocity.  This  is  a  one 
component  system  and  requires  three  scanners.  For  more  information  on  this 
design  see  Econonou  (1986). 


1.3  Objective  of  this  Work 


The  motivation  for  developing  a  rapidly  scanning  LDV  is  to  obtain 
space-time  data  for  and  a  better  understanding  of  separated  turbulent  boundary 
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layers.  The  key  to  this  understanding  is  in  the  instantaneous  behavior  of  the 
large-scale  structures  which  produce  the  large-scale  diffusion  of  turbulence  and 
momentum.  A  three  component  rapidly  scanning  LDV  system  would  provide  the 
necessary  "instantaneous"  velocity  profile  vital  to  understanding  separating  and 
separated  turbulent  boundary  layers. 

The  behavior  of  these  large-scale  structures  observed  through  the 
"instantaneous"  velocity  profiles  would  provide  insight  into  a  number  of 
fundamental  issues  about  separated  turbulent  flows.  Among  these  issues  are 

1.  the  relationship  between  the  large-scale  structures  and  the  on-set  of 
instantaneous  flow  reversal. 

2.  the  relationship  between  velocity  and  pressure  fluctuations. 

3.  the  relationship  between  large-scale  structures  and  the  lag  and  hysteresis  of 
separated  flows. 

The  almost  instantaneous  velocity  profile  provided  by  a  three-component 
scanning  LDV  provides  information  to  determine  the  inflows  and  outflows  of 
turbulent  kinetic  energy  and  momentum  and  the  many  of  the  turbulent 
contributions  to  the  pressure  fluctuations.  Such  separated  flow  data  cannot  be 
obtained  by  hot-wire  arrays  because  flow  reversal  or  by  image-processing  of 
particle  tracking  photographs  because  of  the  three  dimensional  nature  of  such 
flows. 
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The  objectives  of  the  present  work  are  as  follows: 


1.  Develop  a  three-component  rapidly  scanning  laser  Doppler  velocimcter 
(RSLDV)  capable  of  scanning  at  least  40  cm  or  longer  at  high  frequencies. 

2.  Develop  a  data  acquisition  system  for  the  system  above. 


1.4  Outline  of  Thesis 


Chapter  2  deals  with  the  design  of  the  RSLDV.  The  components  of  the 
transmitting  optics  are  explained,  as  well  as  the  problems  with  maintaining  a 
coincident  measuring  volume  over  the  scan  length  and  techniques  investigated  to 
deal  with  this  problem.  The  fringe  patterns  formed  are  also  analyzed  and 
discussed. 

Chapter  3  describes  the  data  acquisition  and  control  hardware  and  logic. 

Chapter  4  gives  the  conclusion  of  this  work  and  makes  recommendations 
for  future  work  and  improvements. 
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2.0  Description  of  RSLDV  Optical  Configuration 


"Instantaneous"  three-component  velocity  profiles  are  needed  in  order  to 
gain  a  better  understanding  of  separated  turbulent  flows.  These  velocity  profiles 
would  provide  insight  into  the  behavior  of  the  large-scale  structures  which  are 
responsible  for  the  large-scale  diffusion  of  turbulence  energy  and  momentum. 

Laser  Doppler  velocimetrv  should  be  used  since  it  provides  low  uncertainty 
measurements  and  does  not  disturb  the  flow  field.  An  LDV  system  to  provide 
the  necessary  "instantaneous"  velocity  profiles  should  be  capable  of  scanning  a 
40  cm  height  at  rates  of  up  to  300  scans/scc.  The  LDV  system  should  provide 
concurrent  measurements  of  all  three  velocity  components  using  an  on-axis  beam 
arrangement,  an  on-axis  system  requires  no  modifications  to  the  facilities  at 
VP1&SU.  The  LDV  system  should  also  allow  for  variations  in  scan  height,  scan 
rate,  location  of  the  measurement  volume,  and  the  fringe  spacing  to  be  made 
easily. 

The  LDV  system  described  below  should  fulfill  these  requirements. 
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2  A  General  Description  of  Setup 


As  shown  in  Figure  1,  the  incoming  beam  from  an  Argon  ion  laser  (/  = 
514.5  nm)  is  directed  by  mirror  (Ml)  to  pass  through  a  20X  laser  beam  expander 
(LBE).  The  expanded  beam  is  then  directed  by  mirror  (M2)  to  pass  through  the 
focusing  lens  combination  (LI  and  L2).  A  dual  water  Bragg  cell  (BC)  splits  and 
frequency  shifts  the  incoming  beam.  In  the  horizontal  plane  an  unshifted  and 
two  shifted  beams,  a  -30  MHz  and  a  +30  MHz,  are  produced.  A  +15  MHz 
beam  is  produced  for  the  vertical  plane.  Mirrors  (M3  and  M4)  provide  a  longer 
path  length  before  mirrors  (M5a,b.c,d)to  allow  easy  separation  of  the  multiple 
beams.  Mirrors  (M5a,b,c,d  and  M6a,b.c,d)  are  used  to  steer  the  beams  at  the 
necessary  angles  onto  the  scanning  mirror  (SM)  and  also  provide  path  length 
equalization  between  the  four  beams.  The  beams  encounter  additional  optics 
after  the  scanning  mirror  (SM)  necessary  to  produce  the  probe  volume.  Figure 
2  shows  the  top  and  side  view  of  the  four  beams  intersecting  to  form  the 
measurement  volume.  Figure  3  shows  the  on  axis  view  of  the  four  beams  forming 
the  measurement  volume.  These  optics  and  the  fringe  patterns  formed  are 
described  in  the  following  sections. 

Using  the  same  optical  setup,  the  laser  can  be  used  in  a  multi-mode 
operation  to  allow  a  two  color  (2  =  514.5  nm  and  488  nm)  LDV  system  with  the 
addition  of  some  color  filters.  Beam  A  would  consist  of  both  wavelengths,  beams 
B  and  D  of  wavelength  X  =  514.5  nm,  and  beam  C  of  wavelength  /  =  488  nm. 
The  on  axis  view  of  the  four  beams  forming  the  measurement  volume  for  this 
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system  is  shown  in  Figure  4.  The  advantage  of  this  arrangment  is  that  there  is 
more  laser  power  in  the  measurement  volume  and  only  four  fringe  patterns  are 
formed.  With  fewer  fringe  patterns  the  Doppler  frequencies  can  vary  more 
widely  in  frequency  without  overlapping. 

2.1.1  Beam  Expansion  and  Focusing  Optics 

The  two  lens  combinations  are  used  in  order  for  the  beams  to  cross  at  their 
waists,  and  to  form  an  appropriate  sized  probe  volume  at  the  desired  location  in 
space. 

The  laser  beam  expander  is  positioned  before  the  focusing  lenses  and  is  an 
Oriel  laser  beam  expander  modified  to  accept  microscope  objectives  as  the  input 
lens,  different  microscope  objectives  can  be  used  to  achieve  different  beam 
expansion  ratios.  The  beam  expander  allows  adjustment  of  the  spacing  betweens 
its  lenses,  thus  determining  the  size  and  divergence  of  the  laser  beam  at  the 
focusing  optics.  This  determines  the  size  of  the  probe  volume,  which  influences 
fringe  visibility  and  laser  power  density  in  the  probe  volume. 

The  focusing  lenses  determine  the  location  of  the  beam  waists  and  consist 
of  positive  achromatic  lenses.  The  input  beam  size  and  divergence  to  this  set  of 
optics  determines  the  size  of  the  probe  volume.  The  focal  lengths  and  spacing 
between  lenses  are  shown  in  Table  1. 

The  design  of  these  two  lens  combinations  depends  upon  several  factors: 
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•  The  size  of  the  probe  volume  necessary  to  maintain  a  coincident  measuring 
volume  (see  2.3  for  more  information). 

•  The  minimum  distance  from  the  last  lens  set  (L4)  to  the  probe  volume. 

Since  Guassian  laser  beams  transform  in  a  manner  different  form  optical 
rays  the  normal  thin  lens  equations  (geometric  optics)  could  not  be  used  without 
large  error  (Yariv,  1985).  Calculation  of  Guassian  beam  propagation  through 
optical  devices  usually  involves  matrix  transformations  or  2-D  Fourier 
transformations,  which  can  be  difficult  to  use.  A  more  straight  forward  but  less 
rigorous  approach  is  given  in  Optics  Guide  4  (1988). 

In  terms  of  input  beam  parameters: 


In  terms  of  output  beam  parameters: 
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[2.4] 


/ 


where 

s  =  object  distance  from  lens 

s"  =  image  distance  from  lens 

/  =  focal  length  of  lens 

ZR  =  Rayleigh  range  of  input  beam 

Z'\  =  Rayleigh  range  of  output  beam 

m  =  magnification  factor 

coq  =  input  beam  waist 

co"0  =  output  beam  waist 

Lens  combinations  arc  handled  by  cascading  the  above  equations;  the  image  from 
the  preceding  lens  is  the  object  for  the  next  lens  in  the  optical  system. 

The  design  of  the  focusing  lenses  is  considered  first  (see  Figure  5).  With 
the  probe  waist  (o)"04)  and  the  location  (</v)  of  the  probe  volume  known  and 
assuming  reasonable  values  for /3  andj^  ,  the  following  equations  result 


*w03 


[2.5] 
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Equations  [2.7]  and  [2.8]  are  iteratively  solved  for  the  spacing  between  the  two 
lenses  (du)  and  the  beam  input  waist  (cu03). 

The  beam  expander  optics  is  based  on  a  reversed  Keplerian  telescope  (see 
Figure  5).  The  equations  for  this  lens  system  reduces  to:  • 


^oi  _  f\ 
w02  fi 


d\2  ~f\  +  /2» 


d2  3  =/2 


[2.9] 


Since  cu"02  =  w0i  and  co0,  is  the  beam  waist  from  the  laser,  selecting/]  determines 
the  value  of/^.  The  beam  profile  through  the  optical  system  shown  in  Figure  5 
is  for  the  l/e2  beam  intensity  location,  which  is  a  hyperbolic  function  in  the 
traverse  direction  (Yariv,  1985). 
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This  optical  configuration  of  expanding  optics  and  focusing  optics  has 
advantages  in  the  flexibility  of  adjusting  the  size  and  location  of  the  probe  volume 
independently. 


2.1.2  Scanner 

The  G325DT  scanner  from  General  Scanning  Inc.  is  a  moving  iron 
galvanometer  with  a  position  transducer  designed  specifically  for  closed-loop 
operation.  This  transducer  operates  by  detection  of  capacitance  variation 
between  the  rotating  armature  and  a  set  of  stationary  electrodes.  The  controller 
includes  a  heater  control  regulating  the  temperature  of  the  scanner.  The  mirror 
is  a  front  surface  mirror  with  DuAg  coating,  and  is  flat  to  one  tenth  of  a 
wavelength  over  the  useable  aperture.  The  shaft  wobble  is  typically  below  5 
arc-seconds  and  the  signal  response  time  is  10  ns. 


2.1.3  Receiving  Optics 

The  receiving  optics  consist  of  two  cylindrical  lenses  (CL1,CL2),  a  long 
front  surface  mirror  (Ml 3),  and  the  PMT,  all  assembled  inside  a  black  box 
having  a  side  door  for  any  adjustments  (Figure  6).  The  first  cylindrical  lens 
(CL1)  gathers  light  and  enlarges  the  image  of  the  measurement  volume  in  the 
streamwise  direction.  The  mirror  deflects  light  onto  the  second  cylindrical  lens 
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(CL2)  which  constructs  the  image  of  the  measurement  volume  on  the  PMT 
aperture  for  any  positions  of  scan.  Both  lenses  and  mirror  can  be  moved 
independently  from  each  other  in  the  vertical  and  horizontal  direction.  For  more 
information  on  the  design  of  the  receiving  optics  see  Antoine  (1985). 

2.2  Horizontal  Plane  Optics 

To  measure  the  velocity  components  in  the  horizontal  plane  (U  and  W), 
three  beams  are  used.  The  optical  configuration  follows  the  same  concept  used 
in  the  designs  of  Chehroudi  and  Simpson  (1983)  and  Antoine  (1985).  Beams  A 
(unshifted),  B  (-30  MHz)  and  D  (  +  30  MHz)  form  the  probe  volume  to  measure 
U  and  W  velocity  components  (see  Figures  1.2,3  and  4). 

Beam  A  is  reflected  from  the  scanner  to  mirrors  (M7,M8.M9a)  to  the 
probe  volume.  Mirrors  (M7,M8)  serve  to  equalize  optical  path  length  from  the 
scanner  to  the  probe  volume.  This  assures  coherence  between  the  three  beams 
and  that  the  beams  scan  the  same  position.  Beam  B  is  reflected  from  the  scanner 
to  mirrors  (M9b,M10,Ml  1)  and  beam  D  to  mirrors  (M9d,M12.M13).  The 
position  and  facing  of  mirrors  (M10,MI  1,M12,M13)  can  be  adjusted  to  allow 
changes  in  the  horizontal  fringe  spacing  and  the  position  of  the  probe  volume 
independently.  Mirrors  (M7.M8)  are  adjusted  so  that  the  path  length  from  the 
scanner  to  the  probe  volume  is  equal  for  all  three  beams. 
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This  optical  arrangement  allows  changes  in  fringe  spacing  and  the  probe 
volume  position  to  be  independently  made.  Also,  since  each  beam  has  its  own  set 
of  optics,  system  alignment  can  be  made  individually  for  each  beam  uncoupled 
from  the  other  beams.  The  path  scanned  by  the  three  beams  forming  the 
horizontal  plane  probe  volume  is  a  vertical  line. 


2.3  Vertical  Plane  Optics 


The  vertical  plane  optics,  necessary  to  measure  the  V  component  of 
velocity,  proved  to  be  the  most  difficult  aspect  in  designing  a  three  component 
rapidly  scanning  LDV  system.  There  were  two  major  problems  that  needed  to 
be  solved: 

•  The  angle  of  intersection  between  the  two  beams  forming  the  vertical  plane 
probe  volume  should  have  a  fringe  spacing  on  the  order  of  a  few  microns. 

•  The  probe  volume  tends  to  scan  an  arc;  coincidence  of  the  horizontal  and 
vertical  probe  volumes  is  difficult  to  achieve. 

The  length  of  the  vertical  plane  probe  volume  ( Lp )  measured  as  the 
projection  along  the  central  beam  (.4)  can  be  found  as  follows 
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[2.10] 


0J!)p  OJ0p 

- L_  - t— 

sin  4>  tan  4> 


where 

4>  =  angle  of  intersection  between  beams  A  and  C 
a>0p  =  the  beam  waist  at  the  measurement  volume 
Coincidence  of  the  probe  volumes  is  maintained  when  the  path  of  the  vertical 
plane  probe  volume  is  located  within  this  area.  Outlined  below  are  several 
methods  investigated  to  solve  the  coincidence  problem. 


2.3.1  Plane  Mirror  Pair 

The  simplest  technique  (and  the  basis  for  some  of  the  other  techniques)  to 
form  the  vertical  plane  probe  volume  is  to  use  two  plane  mirrors  arranged  as 
shown  in  Figure  7.  Beam  C  ( +  15  MHz)  is  reflected  from  the  scanner  at  an  angle  • 

of  60  to  beam  A  to  mirror  (M14),  located  at  a  distance  of  dt  from  the  scanner  and 
angle  0,  to  the  horizontal  plane.  Mirror  (Ml 5),  located  at  d2  from  the  scanner 
and  an  angle  of  02  to  the  horizontal  plane,  reflects  beam  C  to  intersect  beam  A  ® 

to  form  the  probe  volume.  Deriving  ray  equations  for  beams  A  and  C 
For  beam  A: 

• 

y  =  (.v  +  x0)  tan  d  [2-11] 


For  beam  C: 
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at  =  0O  +  0 

[2.12] 

tan  0n  —  tan  0, 
r  V/1  tan  a  —  tan  0, 

[2.13] 

y m\  =  xm  tan  a 

[2.14] 

a  =  20,  -  0O  -  0 

[2.15] 

V  =  20,  -  00 

[2.16] 

jcu1  tan  P  +  <Y,  tan  0O  —  (</,  —  d2)  tan  y  —  d2  tan  02  —  y  w, 

[2.17] 

tan  p  -  tan  02 

y M2  ~  i-xM2  ~  xwi)  tan  P  +y  W1 

[2.18] 

Qr  =  -  (20,  -  20,  -  0O  -  0) 

[2.19] 

T  =  (*  ~  xm)  tan  dr  +  yVf2 

[2.20] 

Equation  [2.1 1]  gives  the  location  of  beam  A  at  any  scan  angle  6  ,  where 
x0  is  the  difference  between  the  path  length  of  beam  A  from  the  scanner  to  the 
probe  volume  and  the  probe  volume  distance  from  the  scanner.  Equations 
[2.13,2.14]  give  the  location  of  beam  C  on  mirror  (M14)  and  [2.17,2.18]  on  mirror 
( M 15).  The  location  of  the  beam  intersection  can  be  found  by  setting  equation 
[2.1 1]  and  [2.20]  equal  to  each  other  and  solving  for  x. 
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[2.21] 


y M2  -  xm  tan  gr  -  x0  tan  0 
tan  0  —  tan  0r 

Equations  [2.21]  and  [2.11]  give  the  .r,y  position  of  the  vertical  probe  volume  as 
the  beams  are  scanned.  From  these  two  equations  a  solution  for  the  Five  variables 
(0O,  0,,  d2,  du  d2)  needs  to  be  found  that  maximizes  the  vertical  and  horizontal 
probe  volume  coincidence.  Since  the  horizontal  probe  volume  scans  a  vertical 
line,  a  solution  can  be  found  by  considering  the  minimization  of  the  following 
integral 


which  can  be  viewed  as  the  difference  squared  between  the  areas  swept  by  the 
two  probe  volumes.  An  optimization  algorithm.  Sequential  Linear  Least  Squares 
Programming  (SLLSQP)  (Kraft,  et  al.  1981),  was  used  to  find  the  minimum  of 
equation  [2.22].  No  global  minimum  was  found  using  this  approach,  even  when 
0O  was  held  constant  at  a  suitable  value,  due  to  the  many  peaks  and  valleys  of 
many  orders  of  magnitude  across  the  domain  space. 

A  brute  force  method  was  then  used  to  find  the  minimum  of  eq.  [2.22];  to 
reduce  computer  time  and  increase  accuracy,  the  above  equations  were 
reformulated  to  reduce  the  number  of  free  variables.  To  eliminate  two  of  the 
variables,  the  beams  were  required  to  have  an  angle  of  intersection  of  4>  and  to 
cross  at  the  required  probe  location  xp  at  a  particular  scan  angle  0.  in  this  case 
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This  enables  02  and  d2  to  be  calculated  in  terms  of  the  other  variables  as 

follows 


02  =  T  «■  -  2e<  +  «o) 


[2.23] 


x2  tan  d2  —  >s  —  *1  tan(20t  —  0O)  -f  yx 
tan  d2  —  tan(20!  -  0O) 


[2.24] 


where 

*1  =  4 

y,  =  tan  d0 

xp  tan  dr{d)  -yp  -  .vv,i(0)  tan  p(0)  +yu,(0) 

tan —  tan  p(6) 

y2  =  (x2  -  xp)  tan  0r{6)  +  yp 

yP  =  {xP  +  x0)  tan  0 

The  minimization  algorithm  increments  d,  across  its  domain  while  holding  0,  and 
90  constant;  the  computer  program  for  this  algorithm  is  in  appendix  A.l.  The 
values  of  d2  and  d2  are  calculated  at  each  d  step  via  cqns.  [2.23]  and  [2.24].  The 
value  of  F  is  calculated  from  eqns.  [2. 11  -  2.22].  The  integral  is  approximated  by 
using  Simpson's  rule.  A  minimum  of  F  over  the  d{  domain  is  obtained  at  each 
0,  step.  The  program  can  be  run  again  with  different  60  and  <£  values,  some  of 
the  results  are  shown  in  Table  2. 

The  advantages  of  this  optical  arrangement  are  that  it  is  simple  and 
flexible:  different  probe  volume  locations  and  fringe  spacing  can  be 
accommodated.  The  disadvantage  is  that  coincidence  is  not  maintained  over  a 
large  scan  length  as  can  be  seen  from  Figure  8. 
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2.3.2  Concave-Plane  Mirror  Pair 


The  path  of  the  vertical  probe  volume  formed  by  the  two  plane  mirrors 
suggests  that  a  concave  mirror  substituted  for  one  of  the  plane  mirrors  would 
tend  to  converge  the  scan  path  closer  to  a  vertical  line.  This  optical  arrangement 
can  be  seen  in  Figure  9.  The  development  of  the  equations  are  similar  to  the 
above  system 


xc  =  dx  —  R  sin  0| 

[2.25] 

T,  =  dx  tan (0„  +  0max)  +  R  cos  0, 

[2.26] 

—  b  +  J  b1  -  4 ac 

X  '«  -  2  a 

[2.27] 

y\ n  =  x\n  tanl',o  + 

[2.28] 

where 

a  -  1  4-  tan:(0o  4-  0) 
b  =  -  2xc  -  2yc  tan(0o  +  d) 
c  =  x2  +  yc2  -  R2 

[2.29] 

[2.30] 


a  =  dn  -t-  Q 


0  “  max 


yj  (d{  —  x  Uj)  4-  (c/j  tan  a  —  \/i )" 
<5  =  2  sin  - — - 
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P  =  20,  ~  0O  -  0ma* 


[2.31] 


7  =  2(0,  -  5)  -  0O  -  e 


[2.32] 


*1/2 


(< d2  —  d{)  tan  P  +  d\  tan  a  —  d2  tan  02  +  xui  tan  y  -  Yxn 

tan  y  -  tan  d2  [--33] 


y M2  =  ixM2  ~  *.wi )  tan  y  4-  ym  [2.34] 

dr  =-  [2(0,  -6)-  202-0o-0]  [2.35] 

The  location  of  the  center  of  the  radius  of  curvature  for  the  concave  mirror  is 
given  by  eqns.  [2.25]  and  [2.26].  Beam  C s  location  on  the  concave  mirror  is  given 
by  eqns.  [2.27]  and  [2.28].  Equation  [2.30]  describes  the  relative  change  in  angle 
with  respect  to  0,  of  the  mirrors  surface  for  the  beam  as  it  scans.  The  beam's 
position  on  the  plane  mirror  is  given  by  eqns  [2.33]  and  [2.34]. 

Variables  d2  and  02  are  again  eliminated  by  using  the  method  outlined  in 
sec.  2.3.1.  This  results  in  the  following 

02  =  -y(4>-20,  +  0o)  [2.36] 


d,  = 


yp  -  xp  tan  0r(0max)  ~T l  +  *i  tan  y(0max) 
tan  y  -  tan  0r(0max) 


[2.37] 


where 

*i  =  d{ 
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y.  =  d\  tan (0O  +  0mJ 


y,  =  (*,  +  Jfo)  0max 

The  beam  intersection  is  found  from  eqn.  [2.21]  and  the  values  above.  The  same 
type  algorithm  is  employed  to  minimize  F  from  eqn  [2.22]  except  another  loop  is  ^ 

added  ic  increment  the  radius  of  curvature  for  the  concave  mirror  (see  Appendix 
A. 2).  Some  results  are  shown  in  Table  3  and  Figure  10  shows  the  probe  volume 
path  for  the  minimum  solution.  % 

This  technique  results  in  better  coincidence  than  that  of  the  two  plane 
mirrors,  but  changes  in  either  the  probe  volume  location  or  the  fringe  spacing 
may  need  a  different  concave  mirror  to  obtain  the  best  coincidence.  ^ 


2.3.3  Prism 


A  prism,  as  shown  in  Figure  11,  was  another  technique  investigated  to 
form  the  vertical  plane  probe  volume.  The  following  equations  were  derived 
considering  the  surface  of  the  prism  as  a  dielectric  interface. 

a  =  d0  +  d  [2.38] 


$1  -  90  -  dp  +  —  +  0O 


+  d 


[2.39] 


d  =  8{  +  sin  '[v'V2  -  sin20,)  sin  e  -  sin  0,  cos  c]  -  e  [2.30] 


d2  =  sin  '(  sin  ) 


[2.41] 
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P  =  ep  -  90  -  j-  +  d2 


[2.42] 


a  e 

[2.43] 

*  =  bp  +  f 

[2.44] 

6r  =  6  +  &g  —  (5 

[2.45] 

D  tan(0o  +  dmax)  -  D  tan  y 
'ril  —  tan  x  —  tan  y 

[2.46] 

>, i  =  xs\  tan  a 

[2.47] 

D  tan(fln  +  0max)  —  ys\  +  xs\  tan  P  -  D  tan  ip 
tan  p  —  tan  t p 

[2.48] 

yS2=ysi  +  ixs2  -  -Yii) tan  « 

[2.49] 

>52  “  -ro  tan  6  -  xs2  tan  dr 
tan  6  —  tan  dr 

[2.50] 

The  same  type  of  algorithm  is  used  to  minimize  F,  sec  appendix  A. 3.  Some 
results  are  shown  in  Table  4.  This  minimization  technique  did  not  work  very  well 
for  this  optical  system,  due  to  oversimplification  in  the  optical  system  parameters. 
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2.3.4  Lens,  Plane  Mirror  Pair 


As  shown  in  Figure  12,  beam  C  passes  through  a  plano-convex  lens  to  the 
plane  mirror  system  described  in  sec.  2.3.1.  This  method  evoked  from 
observation  of  the  prism  system  above,  a  lens  is  used  like  a  prism  but  provides  a 
varying  degree  of  deflection  as  the  beam  is  scanned.  The  equations  for  this 
system  are 


dL  tan(dn  -f  flmax)  -  dL  tan (90  +  +  #max) 

tan(d0  +  0)  -  tan(90  +  0o  +  0maJ 


y i  =  xL  tan(0o  -I-  d) 


—  •  -1  r  xi-  ~  dj  i 
C~Sm  R  sin(0o  +  0max) 


where 


R  =  the  radius  of  curvature  of  the  lens  =f(n-  1) 


5  =  ~d)  +  sin  'l^/n2  -  s\n2(dmax  -  6)  sin  t 


-  sin(0max  -  6)  cos  c]  -  e 


-X  Wl  — 


xL  tan(0o  4-  5  +  6)  —  yL  +  dx(  tan  60  —  tan  0,) 
tan(0n  +  <)  +  0)  —  tan  (?) 
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[2.52] 


[2.53] 


[2.54] 


[2.55] 


J  v/i  =  -xu\  tan(0„  +  5-1 -  6)  +  vL  -  x,  tan(f?0  +  5  +  0)  [2.56] 


P  =  20,  -  0O  -  5  -  0 


[2.57] 


y  =  20,  -  0O 

9 

( d2  -  dx)  tan  y  +  xul  tan  P  +  dx  tan  0O  -  d2  tan  02  —  v,n 

[2.58] 

X,U2  tan  p  —  tan  02 

• 

y.U 2  =  (XM2  ~  x\fl)  tan  P  +>’.vn 

[2.59] 

dr  =  -  (20,  -  202  -  0Q-5-0) 

[2.60] 

The  variables  iU  and  02  are  found  from  eqns  [2.23]  and  [2.24],  except  they  are 
evaluated  at  0  =  0max.  The  same  algorithm  is  used  to  minimize  F,  eqn  [2.22],  with 
the  focal  length  of  the  lens  being  changed  for  each  run  as  necessary.  Results  arc 

<• 

shown  in  Table  5  and  Figure  13. 

This  technique  maintains  the  best  probe  volume  coincidence  but  requires 
^  more  optical  components,  so  set  up  is  a  little  more  complex.  Changes  in  the  probe 

volume  location  or  fringe  spacing  can  be  made  by  changing  the  position  and 
angles  of  the  optics  and  may  also  require  a  different  focal  length  lens  for  best 
m  results. 


2,3.5  Dual  Scanners 

Another  technique  that  can  be  used  to  maintain  the  coincidence  of  the 
probe  volumes  is  to  use  two  scanner.  The  second  scanner  would  be  located  some 
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height,  ys2,  above  the  primary  scanner  and  would  control  the  scan  angle  of  beam 
C  .  The  position  output  of  the  primary  scanner  would  be  used  to  drive  the 
secondary  scanner  through  a  transfer  function.  The  ideal  transfer  function  is 
derived  by  determining  the  angle  that  beam  C  would  have  to  be  directed  to 
intersect  beam  A  at  the  probe  volume,  given  by  Eq.  2.11.  This  is  developed  as 
follows: 

y  =  x  tan  ds2+ys2 


reqno  [2.61]  Setting  Eq.  2.11  and  2.61  and  solving  for  6sl  gives: 

„  .  -ir  ('r  +  *o)  tan  6  —yS2 

ds2  =  tan  L - y - J 


[2.62] 


Figure  14  shows  the  relationship  between  the  scan  angles  for  the  two  scanners 
with  the  current  setup.  The  position  output  from  the  scanners  is  a  voltage 
linearly  dependent  on  the  scan  angle.  Also,  the  scan  angle  is  linearly  dependent 
on  the  input  signal  to  the  scanner  control.  Thus,  the  transfer  function  shown  in 
Figure  14  applies  to  the  position  output  signal  and  input  signal  for  the  two 
scanners. 

This  arrangement  requires  no  addtional  optics,  besides  the  two  scanners, 
and  is  easy  to  setup  and  align.  Probe  vloume  coincidence  can  be  maitained  to  a 
high  degree,  better  than  the  above  purely  optical  methods,  by  this  method 
dependent  on  the  accuracy  that  the  transfer  function  could  be  implemented  in  an 
electronic  circuit.  However,  scanner  drift  and  non-linearities  between  scanner 
angle  and  postion  output  could  be  a  problem. 
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2.4  Doppler  Frequencies  for  the  System 


The  interference  phenomena  between  the  four  laser  beams  produces  six 
moving  fringe  patterns  for  the  single  color  system,  as  shown  in  Figure  3,  and  four 
moving  fringe  patterns  for  the  two  color  system  (see  Figure  4),  that  permit  the 
measurement  of  the  U,  V,  and  W  components  of  velocity.  These  moving  fringe 
patterns  produce  Doppler  frequencies  around  15  MHz,  two  30  MHz,  and  60 
MHz;  the  single  color  system  as  additional  Doppler  frequencies  around  15  MHz 
and  45  MHz.  The  Doppler  frequencies  can  be  determined  by  taking  the  dot 
product  of  unit  vector  of  the  moving  fringe  pattern  with  the  velocity  vector,  this 
results  in  the  following  equations  (see  Figure  2) 

>rd  if~  60MHz)  =  -  U  [2.61] 


I-abV-I 0MHz)  =  — l/cosy  —  (F—  Vs)  sin  -y  sin  6 
+  IF  sin  y  cos  d 


[2.62] 


XAD  if  —  30MHz)  =  —  U  cos  y  +  (V  —  Vs)  sin  —  sin  6 

,  [2-63] 

-  IF  sin  y  cos  6 


XAC{f-  15MHz)  =  {V-  Vs)  cos(  -y  —  d)—  lFsin(y-0)  [2.64] 
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[2.65] 


abc  if -45MHz)  =  -  U 


sin  a 


+  (V-VS) 


sin(0  —  d)  +  cos  a  sin  d 


+  W 


2  sin  iJ/2 
cos(0  —  d)  —  cos  a  cos  Q 
2  sin  c/2 


2  sin  <J/2 


;.CD(/- i5MHz)  =  -<y 


sin  a 


sin(0  —  d)  +  cos  a  sin  d 


-  W 


2  sin  £/2 
cos(0  -  d)  —  cos  a  cos  Q 
2  sin  <5/2 


2  sin  c/2 


[2.66] 


where 


£  =  cos~'(  cos  a  cos  0) 


2  sin  a 


J-ab  —  S-ad  — 

}'AC  — 


AB  AD  2  sin  a/2 


2  sin  0/2 


Ad, - ).rn  — 


'*c  'co  2  sin  m 
The  above  equations  can  be  solved  for  velocity  components  U,  V,  and  W 

in  two  ways.  If  the  frequency  difference  between  equations  2.62  and  2.63  can  be 

obtained  the  following  system  of  equations  result 


V 

v-vs 

= 

w 

-  1  0 

Q  sin(0/2  -  0) 

2  sin  a/2  cos  0/2 

Q  cos(0/2  -  6) 

2  sin  a/2  cos  0/2 


0 


IbdV-  60 MHz) 


cos  6 
cos  0/2 


^abWab  -  /ad) 


[2-67] 


sin  0 
cos  0/2 


;v!C(/-  153/Z/r) 
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The  frequency  difference  can  be  obtained  by  bandpass  filtering  around  30  MHz, 
then  squaring  the  resulting  signal  and  low  pass  filtering,  this  can  be  performed 
analogically  before  sampling  or  digitally  during  processing.  If  the  frequency 
difference  between  the  two  30  MHz  signals  can  not  be  obtained,  then  the 
following  system  of  equations  can  be  solved  for  the  velocity  components 


V 


V  —  Vr 


IV 


0 


cos^_  4-  sin(<ft/2  -  0) 
2  sin  a/2  cos  </>/ 2 


cos(</>/2  —  d) 


cos  0 
cos  </>/ 2 

sin  0 


sin  a/2  cos  </>/ 2  cos  </>/ 2 


60 MHz) 


30 MHz) 


?.AC(f-  15  MHz) 


[2.68] 


In  the  equation  above  the  minus  sign  (-)  refers  to  the  Doppler  frequency  from  eqn 
[2.62]  and  the  plus  sign  (  +  )  to  eqn  [2.63].  For  the  single  color  system  additional 
signals  around  15  MHz  and  45  MHz  (eqns.  2.65,2.66)  are  present.  With  these 
two  equations  and  the  U  velocity  component  calculated  from  one  of  the  system 
of  equations  above,  a  second  set  of  estimates  for  the  V  and  W  velocity 
components  can  be  obtained.  The  two  sets  of  estimates  can  then  be  averaged. 


2.4.1  Uncertainty  Analysis 

An  uncertainty  analysis  can  be  performed  on  the  velocity  estimates 
calculated  from  the  equations  above  using  a  method  given  in  Holman  (1984). 
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We  wish  to  estimate  the  uncertainty  in  the  result  R  on  the  basis  of  the 
uncertainties  in  the  primary  measurements.  The  result  R  is  a  giver  function  of 
independent  variables  x2, ...,  x„ .  Thus, 

R  =  R{xl,x2,  [2.69] 


The  uncertainty  in  R,  wR,  given  the  uncertainty  in  the  independent  variables, 
w,,  vv:, ...,  wn,  is  given  by 


wp  = 


8R  ...  x2  ,  ,  8R  ,2 

Tx rU|) 


+ 


,  (  oR  ,2 


l'/2 


CX„ 


[2.70] 


It  should  be  noted  that  the  uncertainty  in  the  independent  variables  for  this 
uncertainty  analysis  method  are  at  the  same  given  odds  and  that  the  variables 
are  assumed  to  be  normally  distributed.  The  partial  derivatives  in  the  equation  • 

[2.70]  above  can  be  approximated  by 

SR  __  R{x\,x2,...,xi  +  Axj,...x„)-R{xl,x2,...,xi,...x„)  # 

ox,  A.v,  L“'  J 


This  is  useful  when  the  data  reduction  is  rather  complicated  and  the  analytical 
derivatives  are  difficult  to  derive. 

The  uncertainty  analysis  described  above  will  be  applied  to  the  two  color 
RSLDV  setup  with  velocity  estimates  calculated  from  equation  [2.67].  For 
illustrative  purposes  the  following  setup  parameters  will  be  used. 

Beam  intersection  angles: 
a  =  6°  ±  0.5  % 
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<f>  =  10°  ±  0.5  % 

Fringe  spacings: 

).BD  =  2.46  ,um  ±  0. 1  % 

XAB  =  aad  =  4.92  um  ±0.1  % 
aac  -  2.8  /j m  ±0.1  Vo 
Scan  velocity 

Vs  =  60  ±  0.05  m/s 
and  the  scan  angle  will  be  taken  as: 

6  =  0°  ±0.01° 

We  will  consider  a  case  where  the  signal-to-noise  ratio  (SNR)  of  the  PMT  signal 
is  20  dB  and  signal  processing  is  performed  via  the  fast  Fourier  trasnform  (FFT) 
with  zero-padding  and  log-parabolic  interpolation  algorithm  (Shinpaugh,  1989) 
with  1024  points  in  the  data  set.  The  uncertainty  in  the  frequency  estimates 
based  on  this  algorithm  and  SNR  is  ±10  KHz.  The  laser  wavelengths  and  the 
Bragg  cell  frequencies  are  known  with  enough  accuracy  that  they  make  negligible 
contribution  to  the  overall  velocity  uncertainties. 

For  this  analysis  we  consider  the  flowficld  around  a  wing-body  junction 
flow,  from  which  we  can  determine  a  worst  case  and  best  case  senario.  The  worst 
case  will  occur  for  the  following  velocity  componets: 

U  =  ±30  m,s,  V  =  -2  m/S,  and  W  =  16  m  s 
and  the  best  case  will  be  for: 

U  =  0  m/s.  V  =  2ms,  and  W  =  0  m./s 
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The  partial  derivatives  of  the  velocity  estimates  from  eqn.  [2.67]  were  calculated 
analytically  with  respect  to  the  variables  a,  <p,  d ,  abd,  Aab,  Aac,  f  and  I7 ,  . 
Using  the  values  given  above  for  independent  variables  and  their  uncertainties  in 
eqn.  [2.70]  we  can  calculate  the  worst  and  best  case  uncertainties  in  the  velocity 
component  estimates.  For  the  worst  case  we  obtain  the  following  results: 

U  =  30  ±0.04  m/s 

V  =  -2  ±0.1  m/S 
W  =  16  ±0.47  m/s 

and  at  best  the  uncertainties  will  be 
U  =  0  ±0.02 

V  =  2  ±0.1  m/s 
W  =  0  ±0.47  m/s 

However,  through  careful  calibration  it  may  be  possible  to  reduce  these 
uncertainties  in  the  velocity  estimates  to  some  degree. 
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3.0  Data  Acquisition  and  Control  System 


Ideally,  in  a  scanning  LDV  system  one  needs  to  know  the  exact  position 
of  the  measurement  volume  at  the  time  when  a  "validated"  velocity  measurement 
is  detected.  Since  the  PMT  signal  contains  six  Doppler  frequencies  in  high  noise, 
this  precludes  the  use  of  an  online  signal  processing  system.  Therefore,  the  PMT 
signal  needs  to  be  digitized  and  stored,  along  with  the  position  of  the 
measurement  volume,  for  later  processing.  Figure  15  shows  the  general  layout  for 
such  a  data  acquisition  and  control  system,  each  of  the  elements  are  described 
below. 


3.1  IBM  PC-RT  Computer 


The  computer  features  a  32  bit  RISC  (Reduced  Instruction  Set  Computer) 
technology  CPU,  16  Megabytes  of  fast  memory  and  AT  compatible  data  bus  and 
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card  slots.  The  machine  runs  under  A1X  2.2.1  which  is  a  UNIX  based  operating 
system.  Two  400  MB  DASD  (Direct  Access  Storage  Device)  drives  and  a  940 
megabyte  optical  disk  drive  are  connected  to  the  RT  via  a  SCSI  (Small  Computer 
Systems  Interface)  adapter.  The  two  DASD  drives  provide  online  storage  of 
LDV  data  during  data  acquisition  and  the  optical  disk  drive  provides  permanent 
storage  of  the  LDV  data.  Sustained  data  write  rates  of  up  to  410  KB  sec  are 
possible  using  the  DASD  drives  and  the  SCSI  interface. 

A  Jorway  corp.  SCSI  CAMAC  controller  interfaces  the  RT  and  the  two 
transient  recorders.  Data  transfer  rates  of  up  to  1.5  MB, sec  are  possible  over  this 
interface. 


3.2  Transient  Recorders 


Both  transient  recorders,  from  the  LeCroy  corporation,  have  the  TR8828c 
digitizer.  This  digitizer  samples  at  200  mega-samplcs/scc  with  8  bit  resolution 
and  100  MHz  bandwidth.  The  transient  recorder  used  to  record  the  PMT  signal 
has  16  MB  of  memory,  the  other  transient  recorder  has  32  KB  of  memory  and  is 
used  to  record  the  position  output  from  the  scanner  control. 
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3.3  Array  Processor 


The  array  processor  is  a  ZIP  3232-20  from  Mercury  Computer  Systems. 
This  system  features  20  MFLOP  performance,  2  MB  of  fast  memory  and  a  40 
MB/sec  data  way.  Built  in  software  routines  are  callable  from  fortran  and  new 
algorithms  can  be  programmed  via  a  software  interface.  A  1024  point  real 
(non-complex  data)  FFT  (fast  Fourier  transform)  can  be  computed  on  this 
system  in  1.3  msecs. 

3.4  Data  Acquisition  Procedure 

The  whole  data  acquisition  procedure  is  designed  to  minimize  the  amount 
of  time  required  to  gather  a  statistically  significant  record  length,  about  5  seconds 
or  more  of  data.  This  is  done  so  that  variations  in  the  tunnel  properties  can  be 
minimized  during  the  data  acquisition  procedure.  A  total  of  80  transient  recorder 
records  are  taken,  this  is  about  1.25  gigabytes  of  data. 

The  RT  computer  controls  the  entire  data  acquisition  and  transfer  process. 
Both  transient  recorders  are  triggered  at  the  same  time  when  the  measurement 
volume  is  beginning  a  scan.  Data  is  taken  continuously  until  the  transient 
recorder  sampling  the  PMT  signal  has  filled  the  16  MB  of  memory,  which  sends 
a  memory  full  signal  to  the  RT  over  the  GPIB  interface.  The  RT  then  reads  the 
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memory  from  the  transient  recorder  and  writes  this  data  to  a  file  on  the  DASD 
disks;  this  operation  takes  approximately  one  minute  to  complete.  Once  the 
entire  16  MB  of  memory  has  been  transferred,  the  R.T  reads  the  memory  of  the 
other  transient  recorder  containing  the  scan  position  and  writes  this  data  to  a  file. 
The  whole  process  is  repeated  until  the  desired  number  of  records  has  been 
obtained. 

The  files  are  stored  with  a  sequential  type  file  name,  so  the  order  of  the 
files  can  be  maintained.  Also  stored  at  the  beginning  of  each  file  is  the  system 
time  so  that  the  exact  time  of  each  velocity  realization  can  be  easily  determined. 
The  data  can  be  stored  permanently  to  the  optical  disk  once  all  the  data  has  been 
collected;  more  data  can  be  taken  at  another  scan  location  or  the  data  can  be 
processed. 

Using  one  of  the  signal  processing  algorithms  given  by  Shinpaugh  (1989), 
the  data  is  read  from  disk  and  transferred  to  the  array  processor  to  determine  the 
velocity  of  any  Doppler  bursts  in  the  record.  During  processing,  the  computer 
keeps  track  of  the  time  and  location  of  each  velocity  measurement. 


3.5  Seeding 


In  order  to  improve  the  signal-to-noise  ratio  of  the  received  Doppler 
signals  mono-dispersed  polylatcx  spheres  will  be  used  for  seeding.  By  using 
mono-dispersed  particles  the  seeding  level  can  be  adjusted  so  that  there  is  on 
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average  only  one  particle  in  the  measurement  volume  at  one  time,  thus  avoiding 
the  noise  generated  by  having  multiple  submicron  particles  in  the  measurement 
volume  typically  with  other  seeding  techniques.  Also  with  monodispersed 
particles  it  is  possible  to  determine  the  path  through  the  measurement  volume  by 
the  relative  intensity  of  the  signal. 

The  polylatex  spheres  are  mixed  with  either  alcohol  or  water  to  get  a 
mixture  that  is  approximately  0.25%  by  weight.  An  aerosol  is  generated  from 
this  mixture  by  a  seeder  similar  to  the  design  of  Echols  and  Young  (????).  The 
aerosol  is  introduced  into  the  boundary  layer  near  the  begining  of  the  tunnel:  this 
allows  time  for  the  particles  to  dry  before  reaching  the  measurement  volume.  To 
reduce  particle  clumping  the  percent  solids  of  the  solution  can  be  adjusted. 
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4.0  Conclusions  and  Recommendations  for  Future 
Work 


4.  /  Conclusions 


A  rapidly  scanning  threc-vclocity-componcnt  laser  Doppler  vclocimetcr 
has  been  designed,  to  obtain  near  instantaneous  velocity  profiles  in  complicated  • 

three  dimensional  flows. 

The  optical  configuration  allows  the  size  and  location  of  the  measurement 
volume,  as  well  as  all  the  fringe  spacings  to  be  varied  freely.  The  four  beams  have  ® 

independent  lens  and  mirror  sets  to  avoid  coupling,  making  system  alignment 
easier,  changes  made  for  one  beam  do  not  affect  another.  Path  equalization  is 
maintained.  Techniques  to  maintain  coincidence  between  the  vertical  and 
horizontal  probe  volumes  were  investigated.  The  lens,  plane-mirror  pair 
maintained  the  highest  level  of  coincidence  between  the  two  probe  volumes,  and 
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an  appropriate  lens  has  been  manufactured  by  a  local  optician  for  use  with  the 
RSLDV. 

Six  moving  fringe  patterns  at  (2x)  15,  (2x)  30,  45,  and  60  MHz  are 
produced  by  two  Bragg  cells  for  the  single  color  system.  For  the  two  color 
system,  moving  fringe  patterns  at  15,  (2x)  30,  and  60  MHz  are  produced  by  the 
two  Bragg  cells.  The  signal  around  60  MHz  is  a  direct  measurement  of  the  U 
velocity,  and  the  signal  around  15  MHz  and  the  frequency  difference  between  the 
two  30  MHz  signals  are  used  to  calculate  the  V  and  W  velocity  components.  An 
uncertainty  analysis  applied  to  the  system  of  equations  used  to  calculate  the 
velocity  components  with  a  worst  case  senario,  the  uncertainties  in  U,  V,  and  W 
velocity  estimates  are  ±0.04  m/s,  ±  0.1  m/s,  and  ±0.47  m/'s  respectively. 

A  data  acquisition,  control  and  processing  system  has  also  been  designed 
for  use  with  the  RSLDV.  The  PMT  signal  is  sampled  at  200  MB  sec  by  a 
transient  recorder  and  the  position  of  the  scanner  is  also  recorded  by  a  transient 
recorder.  The  data  from  the  transient  recorders  is  transferred  over  a  GP1B 
interface  to  a  mini-computer  for  storage,  the  data  transfer  rate  is  at  300  KB  sec. 
Online  storage  is  provided  for  up  to  1.25  gigabytes  of  data,  this  translates  into  6 
seconds  worth  of  LDV  data.  Permanent  storage  is  provided  by  an  optical  disk 
drive. 
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4.2  Recommendations  for  Future  Work 


Now  that  an  appropriate  design  has  been  made  for  a  rapidly  scanning 
three-velocity-componcnt  laser  Doppler  velocimcter,  additional  optics  can  be 
acquired  and  the  system  setup.  The  final  setup  may  require  some  changes  in  a 
few  of  the  configuration  parameters,  but  these  adjustments  should  be  minor.  The 
quality  of  the  fringe  patterns  can  then  be  observed  and  the  system  can  be 
calibrated  by  using  a  rotating  disk.  The  RSLDV  will  then  be  used  to  investigate 
some  well  documented  turbulent  flows  in  the  Virginia  Tech  Boundary  Layer 
Wind  Tunnel,  including  some  steady  and  unsteady  separated  flows  and  a 
wing-body  junction  flow. 

Improvements  to  the  optical  system  should  be  made  on  increasing  the 
coincidence  between  the  horizontal  and  vertical  probe  volumes.  This  could  be 
accomplished  by  using  a  glass  wedge  of  the  proper  shape,  instead  of  the  current 
lens,  an  elliptical  or  parabolic  mirror  or  by  using  another  scanner  to  scan  the 
upper  beam  with  a  feedback  control  from  the  original  scanner.  By  improving  the 
coincidence  between  the  two  probe  volumes,  the  size  of  the  measurement  volume 
can  be  decreased,  leading  to  improved  resolution  and  higher  SNR  for  the  signal 
from  the  PMT. 

In  an  extension  to  the  current  design  of  the  RSLDV,  a  pockcl  cell  can  be 
placed  just  after  the  Bragg  cells.  This  permits  two  parallel  scan  lines  that  are 
separated  by  a  small  distance,  AX  (about  6  mm),  in  the  freestream  direction.  The 
pockcl  cell  is  switched  on  and  off,  synchronized  with  the  up-and-down  scanning 
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of  the  measurement  volume.  The  measurement  volume  is  scanned  up  at  one  jc, 
location  while  scanning  down  at  an  other  x2  location.  This  permits  information 
on  the  streamwisc  changes  in  the  flow  to  be  obtained.  Two  separate  PMTs  would 
be  needed  in  this  arrangement. 
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Figure  I.  Top  View  of  Optical  Configuration  of  the  RSLDV  System. 
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Figure  4.  On-axis  View  of  (he  Four  lieams  forming  (he  Measurement  Volume  for  (he  Two  Color  System. 
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Figure  8.  Scan  Path  of  Probe  Volumes:  Plane  Mirror  Pair. 
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Figure  II.  Side  View  of  Vertical  Plane  Oplics: 


.  Schematic  Diagram  of  Data  Acquislian,  Control  and  Processing  System. 


APPENDIX  A 


FIFTH  INTERNATIONAL  SYMPOSIUM  ON 
APPLICATION  of  LASER  techniques  to  fluid  mechanics 

AM) 

WORKSHOP  ON  THE  USE  OE  COMPUTERS  IN  PLOW  MEASUREMENTS 


,',</>  9thl2ili.  1990 

Lisbon.  Portugal 

Signal  Processing  Techniques  for  Low  Signal-to-Noise  Ratio 
Laser  Doppler  Velocimetry  Signals 


K.  A  Slunpaught.  R.  L.  Simpsont.A.  L.  WicksJ.andJ.  L.  Flcmingf 
t  Department  of  Aerospace  and  Ocean  Engineering 
jDepanment  of  Mechanical  Engineering 
Virginia  Polytechnic  Institute  and  State  University 
Blacksburg,  VA  24061 


ABSTRACT 

A  variety  of  methods  have  been  developed  to  obtain 
accurate  frequency  estimates  from  laser  Doppler  velocimetry 
signals.  Rapid  scanning  and  fiber  optic  LDV  systems  require 
methods  for  extracting  accurate  frequency  estimates  with 
computational  efficiency  from  data  with  poor  signal-to-noise  ratios. 
These  methods  typically  fall  into  two  general  categories,  ume 
domain  parametric  techniques  and  frequency  domain  techniques. 
The  frequency  domain  approach  is  initiated  by  transforming  the 
Doppler  bursts  in  the  frequency  domain  using  the  fast  Fourier 
transform  (FFT).  From  this  basic  transformation  a  vaAiy  of 
procedures  have  been  developed  to  optimize  the  frequency 
estimation  accuracy.  The  lime  domain  approaches  are  derived  from 
the  parametric  form  of  a  sinusoid.  The  estimation  of  constants  in 
this  relationship  is  performed  to  satisfy  specific  constraints, 
typically  a  minimization  of  a  variance  expression.  A  comparison 
of  these  techniques  is  presented  using  simulated  signals  and  additive 
white  noise.  The  statistical  bias  and  random  errors  for  each  mchiod 
are  presented  from  200  signal  simulations  at  each  condition. 
Frequency  estimation  via  the  FFT  was  found  to  be  superior  to  the 
lime  domain  techniques  studied. 
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cigenfilter  or  characteristic  polynomial 
eigenvector 

autocorrelation  sequence 
frequency,  Hz 

normalized  frequency  f -flf. 

Doppler  frequency  estimate 
frequency  of  FFT  spectral  bin 
sampling  frequency 
fast  Fourier  transform 
exchange  matrix 

number  of  sample  points  in  data  set 
ith  power  spectral  One  from  PSD 
Pisarenko  harmonic  decomposition 
power  spectrum  density 
autocorrelation  coefficient  for  time  lag  i 
autocorrelation  matrix  of  order  M  +  I 
root-mean-square 
signal-to-noise  ratio 
sampling  period 

modified  data  matrix 

spectral  resolution  for  FFT,  A f  “  I IN&i 

sampling  interval 


INTRODUCTION 


The  increased  complexity  of  measured  flowficlds,  the  higher 
degrees  of  required  statistical  information,  and  the  development  of 
more  sophisticated  User  Doppler  velocimetry  systems  (such  as 
rapidly  scanning  LDV  systems)  require  the  processing  of  low 


signal-to-noise  ratio  (SNR)  signals  to  obtain  the  necessary 
information.  Figure  I  shows  a  simulated  LDV  bunt  at  various 
SNR  levels.  We  would  like  to  be  able  to  obtain  accurate  Doppler 
frequency  estimates  for  real  signals  with  SNR  below  20  dB  down  to 
•  10  dB.  Traditional  LDV  signal  processors,  such  as  counten,  have 
proven  to  be  both  versatile  and  accurate  if  the  signal-to-noisc  ratio 
is  high,  usually  well  above  20  dB.  In  order  to  obtain  accurate 
frequency  estimates  at  low  SNR,  frequency  domain 
(non-paramctric)  or  time  domain  parametric  modeling  processing 
must  be  used. 

Processing  in  the  frequency  domain  offers  considerable 
appeal  to  the  experimentalist  from  the  stand  point  of  speed  of 
compulation  and  familiarity.  The  fast  Fourier  transform  (FFT) 
algorithm  has  been  effectively  applied  for  the  past  25  yean  to 
general  signal  processing  and  has  become  a  fixture  in  the  processing 
of  LDV  dtu.  The  FFT  calculates  a  harmonically  related  set  of 
sinusoidal  components  eoual  to  the  number  of  data  points  in  the 
data  set,  each  with  a  tingle  degree  of  freedom.  Ensemble-avenging 
reduces  the  uncertainty  of  the  spectral  estimates;  this  is  typically  not 
available  or  desirable.  To  improve  the  frequency  resolution  of  the 
FFT,  usually  some  type  of  interpolation  scheme  is  used. 

In  time  domain  parametric  methods,  an  a  priori  assumption 
of  the  deterministic  content  of  the  signal  is  made.  The  nature  of  the 
content  is  assessed  based  on  analysis  or  knowledge  of  the  process 
producing  the  signal,  thereby  increasing  the  available  statistical 
degrees  of  freedom.  Parametric  techniques  make  up  two  classes- 
autoregressive  (AR)  based  methods  and  direct  paramocr  estimation 
techniques.  Autoregressive  methods  work  best  with  wideband 

in  noise  and  tend  to  produce  better  spectral  shapes  but  are 
not  as  good  at  locating  spectral  peaks  as  compared  to  direct 
parameter  estimation  methods.  Direct  parameter  estimation 
methods  work  best  for  narrowband  signals  in  noise  (Marple  1987). 
LDV  signals  fall  into  this  classification. 

The  fast  Fourier  transform  (FFT)  algorithm  with  various 
interpolation  schemes,  the  Pisarenko  harmonic  decomposition 
(PHD)  algorithm,  and  the  Eigenvector  (EV)  algorithm  are 
invesugated  for  processing  low  SNR  LDV  signals.  The  FFT  is 
studied  here  since  it  is  a  widely  used  technique  for  promsing  LDV 
signals.  We  would  like  to  know  its  limits  and  capabilities  in  more 
detail  and  study  the  performance  of  the  various  interpolation 
schemes.  Also,  the  FFT  can  be  used  as  a  basis  to  compare  other 
signal  processing  techniques.  The  PHD  is  computationally  (aster 
than  the  FFT  and  theoretically  has  the  ‘ultimate  resolution*  for 
sinusoids  in  white  noise  (Kiy  4  Mtrple  1981),  which  makes  it  an 
attractive  method  for  LDV  signal  processing.  The  Eigenvector 
(EV)  method  is  a  relatively  new  technique  for  frequency  ostimatinn 
which  obtains  high  resolution  of  narrowband  signals  in  noise,  even 
with  small  data  sets  (Marple  1987). 

Each  algorithm  was  evaluated  with  respect  (o  bias  (mean) 
frequency  error,  root-mean-square  (RMS)  frequency  error,  and 
computation  time  as  ■  function  of  frequency  (fraction  of  sampling 
frequency),  number  of  sample  points,  and  signal-to-noise  ratio 
(SNR).  Testing  was  performed  with  computer  simulated  Doppler 
signals  with  specified  frequency  (fraction  of  sampling  frequency), 
number  of  data  points,  and  SNR.  These  signals  were  generated  by 
adding  white  noise  with  Gaussian  distributed  amplitude  to  a 
Gaussian  weighted  sinusoid  to  obtain  a  simulated  signal  with  the 


desired  parameters  (see  Figure  I).  Statistical  information  presented 
(bias  error  and  RMS  error)  is  based  on  200  realizations.  For  each 
realization,  a  statistically  independent  white  noise  sequence  was 
added  to  the  simulated  Doppler  burst.  For  the  evaluation  of  the 
SNR,  the  variance  of  the  white  noise  was  fixed  to  generate  the 
specific  SNR.  The  stgnal-io-noisc  ratio  is  given  by 
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mean  square  of  signal 
variance  of  noise 
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In  this  paper  all  frequencies  have  been  normalized  by  the  sampling 
frequency,  and  all  frequency  errors  are  given  as  percentage  of  the 
sampling  frequency. 


FAST  FOURIER  TRANSFORM 


The  fast  Fourier  transform  is  employed  in  Doppler 
frequency  estimation  by  calculating  the  power  spectral  density 
(PSD)  of  a  Doppler  burst.  In  the  frequency  domain  as  displayed  in 
the  PSD.  the  desired  signal  is  easily  separated  from  the  noise.  The 
resolution  obtained  from  the  PSD  is  inversely  proportional  to  the 
total  sample  period;  N  spectral  lines  are  obtained  for  N  input  data 
points  spanning  the  frequency  domain  of  -0.5  to  0  JS  of  the  sampling 
frequency.  A  common  technique  to  improve  the  estimate  of  the 
Doppler  frequency  of  an  LDV  burst  is  to  (it  a  parabola  through  the 
spectral  peak  and  several  adjacent  spectral  lines  (Kalb  &.  Crosswy 
19S3).  The  motivation  for  using  a  parabola  is  that  it  is  a  symmetric 
function  and  simple  to  evaluate,  but  it  can  also  be  shown  to  be  the 
logical  function  to  use.  The  Fourier  transform  of  a  Gaussian 
function  or  window  is  also  a  Gaussian  function  or  window.  Thus, 
the  Founer  transform  of  a  Gaussian  LDV  burst  results  in  a 
Gaussian  function  centered  at  the  Doppler  frequency.  Taking  the 
natural  loganthm  of  the  spectrum  results  in  a  parabola  centered  at 
the  Doppler  frequency.  Another  interpolation  technique  is  to  lit  a 
Gaussian  curve  to  the  spectral  lines  (Hishida  el  at.  1989). 


The  application  of  time  domain  weighting  functions 
(windows)  to  LDV  data  blocks  have  been  discussed  by  several 
authors  (Bachalo  el  al.  1989,  Matovic  &  Tropea  1990).  The 
application  of  specific  windows  to  sinusoidal  content  is  defined  by 
the  convolution  of  the  frequency  domain  representation  of  the 
sinusoid.  Symmetric  windows  such  as  the  Hanning  window (l  Una). 
when  applied  to  white  noise  dau.  reduce  the  variance  of  the 
spectrum  and  thus  has  been  promoted  as  an  asset  (Harris  1978). 
Since  the  signal  content  of  LDV  dau  has  a  Gaussian  amplitude 
weighting  already,  the  application  of  additional  window's  to  the  dau 
complicates  the  convolution  of  the  signal  content  while  reducing  the 
variance  of  the  noise  content.  Asa  result,  this  study  was  performed 
independent  of  windows  in  addition  to  the  Gaussian. 


The  estimate  for  the  Doppler  frequency  is  found  from  the 
zero-slope  intercept  of  a  parabolic  regression  applied  (Eq.  2)  to  the 
Largest  and  two  adjacent  spectral  lines.  The  parabolic  fit  applied  to 
an  LDV  burst  is  illustrated  in  Figure  2. 
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The  Gaussian  fit  is  applied  in  a  similar  manner  and  is  given  by  (Eq. 
2)  (Hishida  ei  aJ.  1989) 
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A  log  parabolic  fit  can  be  made  by  uking  the  natural  logrithm  of 
the  spectral  lanes  in  Eq.  (2);  this  is  equivalent  to  Eq.  (3).  This 
method  is  unsuble  for  SNR  levels  below  0  dO.  Zero-padding  is  also 
another  technique  used  to  enhance  the  accuracy  of  estimating  the 
frequency  of  spectral  peaks  (Kay  &  Marplc  1981)  and  can  be  used 
in  conjunction  with  any  of  the  above  interpolation  methods.  For 
the  cases  studied  here,  only  N  set  of  zeros  was  added  to  the  dau  set 
to  give  a  total  record  length  of  2N. 

Some  results  for  the  fast  Fourier  transform  with  various 
interpolation  schemes  applied  to  simulated  LDV  signals  are  shown 
in  Figures  2-1.  Figures  3  and  4  show  the  RMS  frequency  error  as 
a  function  of  SNR  and  the  number  of  sample  points,  N,  at  a 
normalized  (by  sampling  frequency)  frequency  of/-  0.225  for  the 
parabolic  fjt  and  Gtuttua  Gt  with  and  without  zero-peddiag.  The 
frequency/ -  0.225  does  not  fall  on  any  spectral  line  for  Ok  number 
of  points,  N,  studied  and  was  selected  for  this  reason.  As  expected. 


a  doubling  in  the  number  of  sample  points  results  in  a  factor  of 
approximately  2  reduction  in  the  RMS  emor  for  each  technique. 
Also,  the  methods  show  a  region  where  the  RMS  error  is  relatively 
constant  with  SNR;  in  this  region  the  noise  power  does  not 
significantly  impact  the  power  content  in  each  spectral  line.  At 
some  point  the  RMS  error  gradually  increases  with  decreasing  SNR 
until  some  limiting  value  of  SNR  is  reached  where  the  RMS 
frequency  error  increases  rapidly.  This  is  due  to  the  noise 
contributing  more  and  mote  power  to  each  spectral  line,  until  even 
in  the  frequency  domain  the  signal  is  buried  in  noise.  This  elTcct  is 
more  pronounced  for  the  cases  with  fewer  sample  points.  Since  the 
power  due  to  noise  is  distributed  approximately  equally  between  the 
spectral  lines  and  the  power  content  in  the  signal  is  only  contained 
in  a  few  spectral  lines,  as  the  number  of  spectral  lines  decreases,  the 
energy  content  of  each  spectral  line  due  to  noise  increases,  but  the 
energy  content  due  to  the  signal  remains  constant.  The  Gaussian 
interpolation  provides  about  a  factor  of  3  decrease  in  the  RMS 
frequency  error  over  the  parabolic  fit  at  high  SNR  but,  at  low  SNR 
the  methods  are  equal.  Zero-padding  improves  the  RMS  error  by 
an  order  of  magnitude  at  high  SNR.  but  is  the  same  at  the  SNR 
limit  Note  that  zero-padduig  does  not  change  the  SNR  limit;  thus, 
the  underlying  resolution  of  the  FFT  is  not  increased  by  this 
method. 


Figures  5  and  6  show  the  RMS  frequency  error  normalized 
by  the  FFT  resolution  tf  as  function  of  signal  frequency  relative  to 
a  PSD  spectral  line  or  FFT  bin  (round  \fN\).  The  cases  shown  are 
for  N  “  128  points  and  SNR's  of  50  dB  and  0  dB.  Since  the  RMS 
error  is  normalized  by  Csf  the  variation  based  on  the  number  of  data 
points  is  negligible,  except  when  the  SNR  limit  has  been  reached  for 
that  particular  N.  Note  that  in  both  cases  the  RMS  error  is 
symmetric  about  a  spectral  line.  For  the  case  of  SNR  —  50  dB,  when 
the  signal  lies  at  ±  0.5  and  0  for  methods  without  zero-padding  and 
±  0.5,  ±  0.25  and  0  for  mchtods  using  zero-padding,  the  RMS  error 
is  nearly  zero,  below  10°.  As  the  SNR  decreases,  the  RMS  error 
becomes  less  dependent  on  relative  frequency  and  eventually  is  flat 
across  the  domain,  except  for  the  parabolic  fit  without 
zero-padding.  These  results  agree  with  those  obtained  by  Hishida 
el  cl.  (1989)  for  (he  Gaussian  fit  without  zero-padding. 


Figures  7  and  8  show  the  frequency  bias  error  normalized 
by  the  FFT  resolution  i\f  as  a  function  of  the  signal  frequency 
relative  to  a  PSD  spectral  line  or  FFT  bin.  Again,  since  the  mean  a 

or  bias  error  has  been  normalized  by  Of .  there  will  be  negligible  ” 

dependence  on  the  number  of  dau  point  N  above  the  SNR  limit. 

The  bias  error  is  antisymmetric  about  zero  relative  frequency.  The 
parabolic  fit  appears  to  be  biased  toward  the  closest  spectral  line  Mb 
and  the  Gaussian  fit  to  the  next  closest  spectral  line.  The  bias  error 
is  nearly  zero  for  relative  frequencies  of  ±  0-5  and  0  for  methods 
without  zero-padding  and  ±0.5,  ±0.25  and  0  for  methods  with 
zero-padding.  The  bias  error  docs  not  change  much  with  SNR. 
except  when  the  SNR  limit  has  been  reached  where  the  bias  error 
increases  rapidly  (Shinpaugh  1989).  These  results  agree  with  those  9 

obtained  by  Bachalo  ei  oL  (1989)  for  the  parabolic  and  Gaussian 
fit  without  zero-padding. 


The  computation  time  is  constant  with  SNR  for  each  case, 
and  these  is  no  difference  between  using  the  parabolic  fit  or 
Gaussian  fit  at  far  as  compuutioo  time  is  concerned.  The  use  of 
zero-padding  approximately  doubles  the  computational  time  over 
the  FFT  algorithms  without  zero- padding  (Shinpaugh  1989). 


PISARENKO  HARMONIC  DECOMPOSITION 


The  Pisarenko  harmonic  dccompostioo  (PHD)  is  a 
parametric  estimation  method  that  assumes  that  the  process 
consists  of  sinusoids  in  additive  white  noise  (Pisarenko,  1973).  The 
sinusoidal  frequencies  and  white  noise  power  are  obtained  from  an 
eigenvalue  analysis  of  the  time-delay  autocorrelation  matrix.  The 
noise  power  is  equal  to  the  minimum  eigenvalue  of  the 
autocorrelation  matrix  of  appropriate  order  and  the  sinusoidal 
frequencies  are  computed  from  the  roots  of  the  characteristic 
polynomial  (eigcnfilter)  associated  with  the  minimum  eigenvalue 
(Hayes  A  dements  198$). 

The  tutocorrelation  sequence  (ACS)  contains  Information 
concerning  the  power  and  frequency  of  each  sinusoid.  Note  that 
the  white  noise  exists  only  in  the  zeroth  lag-  The  autocorrelation 
matrix  formed  from  the  ACS  will  have  the  form 
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Ihc  minimum  eigenvalue  of  R*.,  and  associated  eigenvector  aM., 
correspond  to  the  white  noise  variance  ar.d  the  rcllcction 
coefficients.  respectively.  The  eigenvector  a*.,  associated  with  the 
minimum  eigenvalue  forms  the  characteristic  equation 
h 

d(r)-£/«,.i(*+ Os'*  (5) 

1—0 

where  /((;)  IS  also  known  as  the  eigenfillcr  of  w(/i).  Tlic  roots  of  the 
eigenlllier  will  lie  on  the  unit  circle  at  angles  IrtfJ  for  ls,  ^  H, 
and  hence  determine  the  sinusoidal  frequencies.  For  M  real 
sinusoids  in  wlme  noise,  the  ACS  must  be  known  for  lags  of  0  to 
-M:  2M  roots  wall  be  obtained  from  the  eigenfilter  (Marple  1957). 

Some  of  the  results  for  the  PHD  algorithm  are  shown  in 
Figures  911.  Figure  9  shows  the  RMS  frequeney  error  as  a 
function  of  SNR  and  number  of  sample  points.  The  RMS  error  is 
extremely  low  for  high  SNR  but  increases  almost  linearly  at  a  rate 
of  a  magnitude  per  decade  decrease  in  SNR.  The  bias  ernor  is 
approximately  zero  for  SNR  above  5  dB,  and  increases  as  the  SNR 
tails  below  this  point,  as  shown  in  figure  10.  Also,  the  number  of 
sample  points  does  not  appear  to  be  a  major  factor  in  the  bias  error. 
The  RMS  frequency  error  as  a  function  of  normalized  frequency 
ar.d  number  of  sample  points  for  the  PHD  algorithm  for  the  case 
of  SN  R  -  SO  dO  is  shown  in  Figure  1 1.  The  RMS  error  is  lowest 
for  frequencies  near  half  the  Nyquist  frequency  and  increases  for 
frequencies  near  zero  and  the  NyquisL  The  RMS  ernor  is  also 
sy  mmetric  about  the  half  Nyquist  frequency.  At  other  SNR  levels 
the  shape  of  RMS  error  as  a  function  of  /  is  approximately  the 
same  (Shinpaugh  1989).  Aktar  ti  al.  (1985)  obtained  similar  results 
for  pure  sinusoids  in  white  noise  by  theoretical  analysis  and 
simulation. 


EIGENVECTOR  ALGORITHM 


The  Eigenvector  method  (EV)  is  similar  to  (he  PHD 
algorithm  in  thai  an  eigenanalysis  is  applied  to  an  autocorrelation 
like  matrix.  Consider,  for  example,  a  signal  containing  a  single 
sinusoidal  component  with  additive  noise.  The  signal  vector  space 
would  therefore  span  a  2M  dimensional  subspace  where  M 
represents  the  number  of  sinusoidal  signals.  Tht 
luiocorrclation-like  matrix  is  of  rank  p  where  p>2M  thus  p-2M 
determines  the  noise  subspace.  The  2M  principal  eigenvectors  of 
(he  data  matrix  predominantly  span  the  signal  subspace,  and  the 
singular  values  (eigenvalues)  of  these  principal  eigenvectors  tend  to 
be  larger  than  the  noise  subs  pace  singular  values.  This  allows  the 
separation  of  the  eigenvalues  into  mostly  signal  subspace  and  a 
mostly  noise  subspace,  effectively  enhancing  the  SNR.  Functions 
of  the  vectors  in  cither  the  signal  or  noise  subspace  can  be  used  to 
provide  frequency  estimates  (Marple  1987). 


The  Eigenvector  method  performs  an  eigenanalysis  on  a 
modified  data  matrix  in  the  form  of 
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For  the  noise-free  complex  exponential  sequence 
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where 


A,  -  Ak  exp  /da 
a,  -  exp[e, 


The  modified  data  matrix,  T,.  formed  from  x  will  have  rank  M  and 
M  nonzero  eigenvalues.  A  singular  value  decomposition  is 
performed  on  the  autocorrelation  lie  matrix,  TJT,,  and  results  in 
(he  following 


T),  Tpv,  —  (8) 

for  I  SiSAf.  Tht  remaining  p  —  M  eigenvalues  are  zero.  The 
principal  eigenvectors  of  T/T,  are  a  linear  combination  of  the 
columns  of  a  matrix  composed  of  signal  vectors.  The  eigenvectors 
corresponding  to  the  zero  eigenvalues  are  orthogonal  to  the  M 
signal  subspace  eigenvectors. 

In  the  presence  of  noise,  the  properties  of  the  data  matrix 
are  not  quite  the  same.  The  principal  eigenvalues  associated  with 
the  signal  subspace  are  typically  larger  than  the  p  —  M  noise 
subspace  eigenvalues  and  therefore  can  be  separated.  The 
eigenvectors  associated  with  the  noise  subspace  are  theoretically 
orthogonal  to  the  sinusoidal  signal  vectors.  Therefore,  the  inner 
product  with  arbitary  weighting  will  be  zero  when  a  noise  vector 
corresponds  to  a  sinusoid  vector.  Inverting  the  inner  product  drives 
the  coincidence  of  the  noise  vectors  and  signal  vectors  to  infinity, 
thus  allowing  the  identification  of  the  frequency  content.  For  a 
more  derailed  description  of  this  algorithm,  see  Marple  (1987). 

Some  results  for  the  EV  method  are  shown  in  Figures  12-15. 
Due  to  the  complexity  of  this  method,  the  number  of  sample  points 
had  to  be  limited  to  below  512  points.  The  RMS  frequency  error 
as  a  function  of  SNR  aqd  the  number  of  sample  points,  N,  for  a 
normalized  frequency  of/ -0.225  is  shown  in  Figure  12.  The  RMS 
error  is  fairly  constant  to  about  20  dB  and  then  rises  at  an  order  of 
magnitude  per  decade.  The  apparent  better  performance  of  the 
method  with  N— 64  points  is  probably  due  to  the  effects  of 
sampling  density.  The  bias  error  as  a  function  of  SNR  and  N  is 
shown  in  Figure  13.  the  bias  error  is  fairly  constant  down  to  about 
10  dB  SNR  but  increases  rapidly  below  -10  dB  SNR.  The  bias  error 
is  not  identically  zero  because  of  the  use  of  a  FFT  to  produce  a 
spectrum  and  the  subsequent  application  of  a  parabolic  Ct  to 
determine  the  frequency.  The  effects  of  normalized  frequency  on 
the  RMS  error  is  shown  in  Figure  14  for  an  SNR  of  50  dB..As  can 
be  seen  from  this  plot,  the  RMS  error  is  symmetric  about/  —  0.25 
,  and  the  sampling  density  has  a  great  effect  on  the  RMS  error.  The 
EV  method  behaves  in  a  similiar  manner  to  the  PHD  algorithm, 
which  it  not  surprising  since  both  use  an  eigenanalysis  on  an 
autocorrelation-like  matrix. 


CONCLUSIONS 

The  RMS  frequency  error  and  the  bias  frequency  errpr  for 
the  FFT  methods  do  not  depend  on  /  except  on  where  /  falls 
relative  to  a  PSD  spectra1  line.  Also,  the  RMS  frequency  error 
becomes  independent  of  /  when  the  SNR  is  below  10  dB.  With 
cither  the  parabolic  or  Gaussian  interpolation  schemes,  the 
zero- padding  procedure  for  the  FFT  reduces  the  RMS  frequency 
error  to  a  few  percent  at  low  SNR  and  an  order  of  magzjtudc  lower 
at  high  SNR  compared  to  the  tame  interpolation  scheme  without 
zero-padding.  Also,  zero-padding  reduces  the  bias  error  typically 
an  order  of  magnitude.  However,  zero-padding  increases  the 
computation  time  by  a  factor  of  approximately  two.  The  Gaussian 
interpolation  method  for  the  FFT  provides  a  factor  of  3 
improvement  in  the  RMS  and  bias  error  at  SNR  levels  above  0  dB 
compared  to  the  FFT  with  the  parabolic  fit.  The  Gaustian 
interpolation  method  with  zero-padding  gives  the  highest  resolution 
for  u grials  with  an  SNR  above  -6  dB.  Below  this  point, 
zero-padding  does  not  improve  the  results,  so  the  Gaussian  Ct 
without  zero-padding  is  to  be  used. 

The  PHD  algorithm  is  typically  an  order  of  magnitude  faster 
computationally  compared  to  the  FFT  algorithms.  The  bias  error 
for  the  PHD  at  higher  SNR  levels  is  smaller  than  tlx  bias  errors 
obtained  for  (he  FFT  algorithms  without  zero-padding  for 
frequencies  tway  from  zero  and  the  Nyquist  frequency.  Only  at 
high  SNR  levels,  above  30  dB,  and  only  at  certain  frequencies  docs 
the  PHD  algorithm  outperform  the  FFT  algorithms  with  respect  to 
the  level  of  RMS  error.  The  EV  method  does  not  compare  well  to 
either  the  FFT  or  PHD  algorithms  except  maybe  where  the  number 
or  temple  points  are  very  few,  below  N- 32-  This  can  be  teen  fat 
Figure  15,  which  shows  the  domain  or  region  at  which  each 
algorithm  is  the  most  effective.  As  an  example,  if  the  SNR  levels 
from  the  photomultiplier  tube  (after  filtering)  were  at  20  dB  then  the 
FFT  with  zero-padding  and  log  parabolic  interpolation  would 
provide  the  most  accurate  frequency  estimates.  As  can  be  teen 
from  Figure  15,  the  switch  pom  is  between  the  two  FFT  algorithms 
occurs  at  the  tame  SNR  level  regardless  of  the  frequency.  However, 
the  dividing  line  between  when  to  use  the  PHD  algorithm  and  the 
FFT  with  zero-padding  and  Gaussian  interpolation  algorithm  is 
dependent  on  the  frequeney  and  the  number  of  samples  m  the  data 
record.  From  this  map  it  is  possible  to  tell  which  algorithm  to  use 
if  an  approximate  SNR  and  frequency  is  known  for  (he  signal  to 
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be  processed. 
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figure  2.  Illustration  of  the  parabolic  fit  to  a  burst  PSD. 
10* 


Parabolic  fit 


10  20 
SNR  (d9) 


Figure  3.  RMS  frequency  error,  n.  SNR  far  the  penboDc  fit 
with  and  without  i*n>- padding  at/  •  02X5. 
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Figure  I.  Simulated  LDV  bunt  dgnab  at  rarlouz  SNR  keels. 


Figure  4.  RMS  frequency  error. r*.  SNR  for  the  Gaimtaa  St 
with  and  without  zero-padding  at/  “  OJli. 
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Figure  S.  RMS  frequency  error  normalized  by  FFT  spectral 
resolution  as  a  function  of  the  signal  frequency  relative  to  a 
spectral  line  for  the  FFT  methods,  SNR  “  50  dB. 


Figure  6.  RMS  frequency  error  normalized  by  FFT  spectral 
resolution  as  a  function  of  the  signal  frequency  relative  to  a 
spectral  Une  for  the  FFT  methods,  SNR  m  0  dB. 
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Figure  7.  Frequency  bUs  error  normalized  by  FFT  spectral 
resolution  u  a  function  of  the  signal  frequency  relative  to  I 
spectral  Une  for  the  FFT  methods,  SNR  -  SO  dB. 


Fit  ore  (.  Frequency  bias  error  normalized  by  FFT  spectral 
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spectra]  line  for  the  FFT  methods,  SNR  -  0  dB. 
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Figure  10.  Frequency  Was  error  n  SNR  for  PHD  at  /  -  (L2S. 


Figure  9.  RMS  frequency  error  vs.  SNR  for  PHD  at  /  -  0  74, 
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Figure  1 1 .  RMS  frequency  error  vs.  frequency  for  PHD  at  SNR 

-  50  dB. 
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13.  Frequency  bias  error  T»  SNR  for  EV  at  /  -  0.225. 
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Figure  15.  SNR  rt  frequency  Identifying  the  tppropriste 
technique  to  use  relillre  to  SNR  sad  frequency.  :  CF  - 
CsunUn  Fit,  ZP  -  zero-padding. 
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Figure  12.  RMS  frequency  error  tj  SNR  for  EV  at /  “  0L22S. 
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Figure  14.  RMS  frequency  error  rs  frequency  far  EV  st  SNR 
-  50  dB. 


